#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#include "../include/var_struc.h"
#include "../include/finite_volume.h"


//Initialize conserved quantities.
void cons_qty_init(const struct cell_var * cv, const struct flu_var * FV)
{
    const int dim = (int)config[0];
    const int num_cell = (int)config[3];
    const double q_0 = config[120];
    for(int k = 0; k < num_cell; k++)
	{
	    cv->U_rho[k]   = FV->RHO[k];
	    cv->U_gamma[k] = FV->RHO[k] * FV->gamma[k];
		
	    cv->U_e[k]     = FV->P[k]/(FV->gamma[k]-1.0) + FV->RHO[k]*q_0*FV->PHI[k] + 0.5*FV->RHO[k]*FV->U[k]*FV->U[k];
	    cv->U_u[k]     = FV->RHO[k] * FV->U[k];			
	    if (dim > 1)
		{									
		    cv->U_v[k]  = FV->RHO[k] * FV->V[k];
		    cv->U_e[k] += 0.5*FV->RHO[k]*FV->V[k]*FV->V[k];
		}
	    if (dim > 2)
		{									
		    cv->U_w[k]  = FV->RHO[k] * FV->W[k];
		    cv->U_e[k] += 0.5*FV->RHO[k]*FV->W[k]*FV->W[k];
		}
	    if ((int)config[2] == 2)
		{									
		    cv->U_phi[k] = FV->RHO[k] * FV->PHI[k];			
		    cv->U_e_a[k] = FV->Z_a[k] * FV->P[k]/(config[6]-1.0) + cv->U_phi[k]*q_0*FV->PHI[k] + 0.5*cv->U_phi[k]*(FV->U[k]*FV->U[k]+FV->V[k]*FV->V[k]);
		}
	}
}


int cons2prim(struct i_f_var * ifv)
{
    const int dim = (int)config[0];
    const double eps = config[4];
    const double q_0 = config[120];
    double phi_e_a, phi_e_b;
	
    ifv->RHO   = ifv->U_rho;
    ifv->U     = ifv->U_u/ifv->U_rho;
    if (dim > 1)
	{
	    ifv->V  = ifv->U_v/ifv->U_rho;
	    if (isnan(ifv->V) || isinf(ifv->V))
		return 0;
	}
    if ((int)config[2] == 2)
	{
	    ifv->PHI = ifv->U_phi/ifv->U_rho;
	    ifv->Z_a = 1.0;
	    /*
	    phi_e_a  = ifv->U_e_a-0.5*ifv->U_phi*(ifv->U*ifv->U+ifv->V*ifv->V);
	    //phi_e_b  = (ifv->U_e-ifv->U_e_a)-0.5*(ifv->U_rho-ifv->U_phi)*(ifv->U*ifv->U+ifv->V*ifv->V);
	    phi_e_b  = ifv->U_e-0.5*(ifv->U_rho)*(ifv->U*ifv->U+ifv->V*ifv->V)-phi_e_a;
	    ifv->Z_a = phi_e_a*(config[6]-1.0)/(phi_e_a*(config[6]-1.0) + phi_e_b*(config[106]-1.0));
	    */
	    if (isnan(ifv->Z_a)||isnan(ifv->PHI))
		return 0;
		else if (ifv->Z_a<(-eps))
		{
		    printf("\n Z_a=%.10lf,\t phi_a=%.10lf\n",ifv->Z_a,ifv->PHI);
		    //ifv->Z_a = 0.0;
		    //ifv->U_e_a = 0.5*ifv->U_phi*(ifv->U*ifv->U+ifv->V*ifv->V);
		}
	    else if (ifv->Z_a>(1.0+eps))
		{
		    printf("\n Z_a=%.10lf,\t phi_a=%.10lf\n",ifv->Z_a,ifv->PHI);
		    //ifv->Z_a = 1.0;
		    //ifv->U_e_a = ifv->U_e-0.5*(ifv->U_rho-ifv->U_phi)*(ifv->U*ifv->U+ifv->V*ifv->V);
		}
	    else if (ifv->PHI<(-eps)||ifv->PHI>(1.0+eps))
		{
		    printf("\n Z_a=%.10lf,\t phi_a=%.10lf\n",ifv->Z_a,ifv->PHI);
		}
	    // ifv->gamma = (phi_e_a*config[6] + phi_e_b*config[106])/(phi_e_a+phi_e_b);
	    ifv->gamma = 1.0/(ifv->Z_a/(config[6]-1.0)+(1.0-ifv->Z_a)/(config[106]-1.0))+1.0; 
	}
    ifv->P     = (ifv->U_e - ifv->RHO*q_0*ifv->PHI - 0.5*(ifv->U_u*ifv->U_u)/ifv->U_rho) * (ifv->gamma-1.0);	
    if (dim > 1)
	{
	    ifv->P -= (0.5*(ifv->U_v*ifv->U_v)/ifv->U_rho) * (ifv->gamma-1.0);
	}
    if (dim > 2)
	{			
	    ifv->W  = ifv->U_w/ifv->U_rho;
	    ifv->P -= (0.5*(ifv->U_w*ifv->U_w)/ifv->U_rho) * (ifv->gamma-1.0);
	    if (isnan(ifv->W) || isinf(ifv->W))
		return 0;
	}
	
    if (isnan(ifv->RHO + ifv->U + ifv->P) || isinf(ifv->RHO + ifv->U + ifv->P) || ifv->RHO < eps)
	return 0;
    else if (ifv->P < eps)
	{
	    printf("P=%.10lf\n",ifv->P);		
	    ifv->U_e = 0.5*(ifv->U_u*ifv->U_u)/ifv->U_rho + eps/(ifv->gamma-1.0);
	    if (dim > 1)			
		ifv->U_e += 0.5*(ifv->U_v*ifv->U_v)/ifv->U_rho;	
	    //ifv->U_e_a = 0.5*ifv->U_phi*(ifv->U*ifv->U+ifv->V*ifv->V) + ifv->Z_a*eps/(config[6]-1.0);	
	}
	
    return 1;
}


int cons_qty_update(const struct cell_var * cv, const struct mesh_var * mv,
		    const struct flu_var * FV, const double tau)
{
    const int dim = (int)config[0];
    const int num_cell = (int)config[3];
    const double K = config[121];
    const double T_c = config[122];
    struct i_f_var ifv, ifv_next;
    double k_rho, k_rho_fix;
    int ** cp = mv->cell_pt;
	
    int p_p, p_n;
    double length;
    int i;
    for(int k = 0; k < num_cell; ++k)
	{
	    k_rho_fix = 0.0;
	    cons_qty_copy_cv2ifv(&ifv, cv, k);
	    for(int j = 0; j < cp[k][0]; j++)
		{
		    if(j == cp[k][0]-1) 
			{
			    p_p=cp[k][1];
			    p_n=cp[k][j+1];
			}				  
		    else
			{
			    p_p=cp[k][j+2];
			    p_n=cp[k][j+1];
			}
		    if (dim == 1)
			length = cv->n_x[k][j];
		    else if (dim == 2)
			length = sqrt((mv->X[p_p] - mv->X[p_n])*(mv->X[p_p]-mv->X[p_n]) + (mv->Y[p_p] - mv->Y[p_n])*(mv->Y[p_p]-mv->Y[p_n]));
					
		    cv->U_rho[k] += - tau*cv->F_rho[k][j] * length / cv->vol[k];
		    cv->U_e[k]   += - tau*cv->F_e[k][j]   * length / cv->vol[k];	
		    cv->U_u[k]   += - tau*cv->F_u[k][j]   * length / cv->vol[k];
		    if (dim > 1)
			cv->U_v[k] += - tau*cv->F_v[k][j] * length / cv->vol[k];
		    if (dim > 2)
			cv->U_w[k] += - tau*cv->F_w[k][j] * length / cv->vol[k];
		    if ((int)config[2] == 2)
			cv->U_phi[k] += - tau*cv->F_phi[k][j] * length / cv->vol[k];
		    if(!isinf(config[60]))
			cv->U_gamma[k] += - tau*cv->F_gamma[k][j] * length / cv->vol[k];
		    k_rho_fix += 0.25*K*cv->rho_phi[k][j];
		}
	    k_rho = k_rho_fix * isgreater(ifv.P/ifv.RHO,T_c);
	    cv->U_phi[k] += - tau*k_rho;
	    cons_qty_copy_cv2ifv(&ifv_next, cv, k);
	    cv->U_phi[k] +=   tau*k_rho;
	    k_rho = k_rho_fix * isgreater(0.5*(ifv.P/ifv.RHO+ifv_next.P/ifv_next.RHO),T_c);
	    cv->U_phi[k] += - tau*k_rho;
	}
    return 0;
}
